home *** CD-ROM | disk | FTP | other *** search
/ Ian & Stuart's Australian Mac 1993 September / September 93.iso / Archives / Applications / Calculators / NumberCrunch 1.41 / Number Crunch II Demo / Library Routines / Text Examples / Fourier Transform < prev    next >
Encoding:
Text File  |  1992-09-03  |  6.6 KB  |  195 lines  |  [TEXT/NCII]
  1. ##########
  2. # Fourier Transforms
  3. #
  4. #   function FFT(Cx)  #  Returns the discrete fourier transform of Cx.
  5. #   function InvFFT(Ck) # Returns inverse discrete fourier transform.
  6. #
  7. #   function SwapSides(V)  # Returns array V with left and right halfs swapped.
  8. #
  9. # 
  10. #  With i=sqrt(-1), the transform is
  11. #
  12. #  Cx=InvFFT(Cq)    # (Sign= +1 transform) :
  13. #    Cx[q] =  Sum p=0…N-1 { exp( +i 2π p q/N ) Ck[p]  }
  14. #  Cq=FFT(Cx)         # (Sign= -1 transform) :
  15. #    Ck[p] =(1/N) Sum q=0…N-1 { exp( -i 2π p q/N ) Cx[q] }
  16. # 
  17. #  which has the normalization
  18. #   (1/N)  sum( |Cx|^2 ) = sum( |Ck|^2 ).
  19. #
  20. # 
  21. #  This text file explains and gives examples of
  22. #  some of the routines in the file All Library Routines, which
  23. #  should be Opened before trying any of these examples.
  24. #
  25. ##########
  26.  
  27.  
  28. #########
  29. # Here is the xCOD:
  30. xfft
  31.   # xCO2 xFFT(zN, Sign:real; cmplxData:ComplexArray; Error:real) , does a 1D Fast Fourier Transform on cmplxData[1..N]. zN must be a power of 2, else Error will be set to -1.; an external program.
  32.  
  33.  
  34. #################
  35. # Interface files:
  36. #
  37.  
  38.   function FFT(Cx)  #  Returns the discrete fourier transform of Cx.
  39. .   var n, Ck, err
  40. . #   Input:     Cx = real or complex array, 2^N points.
  41. . #   Output:   FFT = complex array = discrete fourier transform of Cx.
  42. .   begin
  43. .     Ck = Cx
  44. .     n = size(Ck)
  45. .     xFFT(n,-1,Ck,Err)
  46. .     if Err<>0 then
  47. .       begin 
  48. .         print(" •• ERROR in xFFT, err="+Err)
  49. .         FFT = 0/0  # which is "Not a number"
  50. .       end
  51. .     else
  52. .     FFT = Ck
  53. .   end
  54. .   
  55.  
  56.   function InvFFT(Ck) # Returns inverse discrete fourier transform.
  57. .   var n, Cx, err
  58. . #   Input:   Ck = real or complex 2^N array
  59. . #   Output:  InvFFT = inverse discrete fourier transform, a complex array.
  60. .   begin
  61. .     Cx = Ck
  62. .     n = size(Cx)
  63. .     xFFT(n,1,Cx,Err)
  64. .     if Err<>0 then
  65. .       begin 
  66. .         print(" •• ERROR in xFFT, err="+Err)
  67. .         invFFT = 0/0  # which is "Not a number"
  68. .       end
  69. .     else
  70. .     InvFFT = Cx
  71. .     end
  72. .     
  73.  
  74. #
  75. # The indeces are cyclic in the discrete fourier transform, with
  76. # Ck[N] the same as Ck[0]. FFT expects the "0 frequency" channel to 
  77. # be the first number in the Ck array, which means that 
  78. # the data is stored in the order
  79. #     Ck[0], Ck[1], Ck[2],... , Ck[N-1]. 
  80. # Because of the periodicity, Ck[N-1] is the same what would be 
  81. # Ck[-1]; therefore, if its often convenient to plot Ck and Cx in the
  82. # order Ck[N/2], Ck[N/2+1],... ,Ck[N-1], Ck[0], Ck[1],.. Ck[N/2 -1]
  83. # which puts 0 in the middle.  
  84.  
  85. # Here's a routine to swap the two sides of an array, so that
  86. # the x=0 or k=0 channel will be in the middle when we plot it.
  87.  
  88. # SwapSides(V) exchanges the left and right halves of a vector V.
  89. # Evaluate this by selecting it (or putting the cursor at the last
  90. # period) and hitting enter.
  91. #
  92.   function SwapSides(V)  # Returns array V with left and right halfs swapped.
  93. .    var n, n0, SS, left,right
  94. . # Input:   V[1…N]
  95. . # Output: SwapSides = { V[N/2], V[N/2+1], …, V[1],V[2,],…,V[N/2-1] }
  96. .   begin
  97. .     n = Size(V);  # the length of V
  98. .     n0 = FirstIndex(V); # V is V[n0...n0+n-1]
  99. .     SS = V  # save space for answer
  100. .     left = n0…(n0+n/2+1)  # left half indeces
  101. .     right = (n0+n/2)…(n0+n-1) # right half indeces
  102. .     SS[left] = V[right]
  103. .     SS[right] = V[left]
  104. .     SwapSides = SS    # return the answer
  105. .   end
  106. .   
  107.  
  108. ###########
  109. #  Example:
  110. #
  111.  
  112. N = 128;            # the number of points in the arrays.
  113. Cx[0…N-1] = 0;     # initialize the array.
  114. x = 0 … N-1;     # initialize the x co-ords of the array.
  115. xSwap = -N/2 … N/2-1;     # x co-ords of the "swapped" array.
  116.  
  117. # Here's a "hatbox" , with C=1 for |x|<=4.
  118. Cx[0] = 1;      
  119. Cx[1…4] = 1;
  120. Cx[N-4 … N-1] = 1;
  121. # To see this, plot [xSwap, SwapSides(Cx)]
  122.  
  123. # Finally, do the transform.
  124. Ck = FFT(Cx)           # <<<<<<<<<<<<   This does the actual transform.
  125.  
  126. # Several things went on there behind the scenes.
  127.  
  128. # First, Err didn't exist when xFFT was called, so it was created
  129. # as a single real number.
  130. Err
  131.   Err = 0 
  132.  
  133. # Second, Ck was a real array when passed to xFFT, but xFFT
  134. # expected a complex array. (See the argument list returned
  135. # when FFT was evaluated.)  So, Ck was converted to a complex
  136. # array (with 0's for the complex part), transformed, and returned.
  137. # Although you usually won't need to tell if an array is real or 
  138. # complex (and in fact NCII will at times convert it back and forth
  139. # like this if it needs to and the the imaginary part is zero), you can
  140. # tell how big an array really is by using Size and SizeAsIfReal, 
  141. # which aren't the same if the array or number is complex.
  142.  
  143. # If you really want to see when something's complex,
  144. # here's a function that will beep when it is..
  145. program BeepIfComplex(z)  # beeps if Z is complex.
  146. .   begin
  147. .      if Size(z) <> SizeAsIfReal(z) then beep
  148. .   end
  149. .   
  150.  
  151. # BeepIfComplex(Cx)   # This won't beep when evaluated without the "#".
  152. # BeepIfComplex(Ck)   #  This will.
  153.  
  154. # To see the transformed array, plot [xSwap, SwapSides(Ck)] .
  155. # This plots the real part of Ck. To see the imaginary part, you could
  156. # plot [xSwap, SwapSides(Imaginary(Ck))] ; 
  157. # but since Cx is symmetric, the imaginary part is zero
  158. # (except for numerical errors), so this is boring in this case.
  159. maximum(| Imaginary(Ck) |)
  160.    2.8649365e-16 
  161.  
  162. # Here is the check of the normalization:
  163. sum( |Ck|^2 )
  164.   0.0703125 
  165. 1/N sum( |Cx|^2 )
  166.   0.0703125 
  167.  
  168. # And finally, here's a table of some of the values. Note that these
  169. # have not been hit with SwapSides, so x=0 is the first number.
  170. # Table formats imaginary quantities as two colums of reals, 
  171. # so the last two columns are both Ck.  
  172. #
  173. # The actual table command is commented so that it won't do
  174. # anything if this whole file is evaluated. Just select it and hit
  175. # enter to make the table.
  176. #   j = 0…12; Table(j, Cx[j],Ck[j])
  177. #
  178. # x                 Cx                 Ck
  179. # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
  180.    0                  1                  0.07031       0
  181.    1                  1                  0.06975       -5.455e-17
  182.    2                  1                  0.06807       -5.3077e-17
  183.    3                  1                  0.06534       -1.0133e-16
  184.    4                  1                  0.06161       -4.7384e-17
  185.    5                  0                  0.05701       -8.6654e-17
  186.    6                  0                  0.05165       -7.7217e-17
  187.    7                  0                  0.04568       -1.0006e-16
  188.    8                  0                  0.03928       -2.7712e-17
  189.    9                  0                  0.0326         -4.366e-17
  190.    10                0                  0.02583       -3.173e-17
  191.    11                0                  0.01913       -2.9912e-17
  192.    12                0                  0.01269       -8.5849e-18
  193.  
  194.  
  195.